By most standards exploration geophysics is a relatively young science, with some of the earliest work in the subject area dating back to the 1920's and the renowned CMP approach dating from only the 1950's. In the years since its genesis, however, it has become the oil industry's preeminent approach to finding subsurface petroleum deposits. Although exploration geophysics generally encompasses the three broad subject areas of gravity, magnetics, and seismic, today it is the seismic method that dominates almost to the point of excluding the other disciplines. In fact, a simple count of the number of seismic crews in the field has become one accepted measure of the health of the entire oil industry.
A seismic survey represents an attempt to map the subsurface of the earth by sending sound energy down into the ground and recording the "echoes" that return from the rock layers below. The source of the down-going sound energy might come from explosions or seismic vibrators on land, and air guns in marine environments. During a seismic survey, the energy source is moved across the land above the geologic structure of interest. Each time the source is detonated, it is recorded at a great many locations on the surface of the earth. Multiple explosion/recording combinations are then combined to create a near continuous profile of the subsurface that can extend for many miles. In a two-dimensional (2-D) seismic survey, the recording locations are generally laid out along a straight line, whereas in a three-dimensional (3-D) survey the recording locations are distributed across the surface in a grid pattern. In simplest terms, a 2-D seismic line can be thought of as giving a cross sectional picture of the earth layers as they exist directly beneath the recording locations. A 3-D survey produces a data "cube" or volume that is, at least conceptually, a 3-D picture of the subsurface that lies beneath the survey area. Note that it is possible to extract individual 2-D line surveys from within a 3-D data volume.
A seismic survey is composed of a very large number of individual seismic recordings or traces. In a typical 2-D survey, there will usually be several tens of thousands of traces, whereas in a 3-D survey the number of individual traces may run into the multiple millions of traces. A seismic trace is a digital recording of the sound energy reflecting back from inhomogeneities in the subsurface, a partial reflection occurring each time there is a change in the acoustic impedance of the subsurface materials. The digital samples are usually acquired at 0.004 second (4 millisecond) intervals, although 2 millisecond and 1 millisecond intervals are also common. Thus, each sample in a seismic trace is associated with a travel time, a two-way travel time in the case of reflected energy. Further, the surface position of every trace in a seismic survey is carefully recorded and is generally made a part of the trace itself (as part of the trace header information). This allows the seismic information contained within the traces to be later correlated with specific subsurface locations, thereby providing a means for posting and contouring seismic data, and attributes extracted therefrom, on a map (i.e., "mapping"). The signal that is sent down into the earth is called a seismic waveform or seismic wavelet. Different seismic waveforms are generated depending on whether the source is an air gun, dynamite or vibrator. The term "source signature" or "source pulse" is generally used to describe the recorded seismic character of a particular seismic waveform.
A seismic source generated at the surface of the earth immediately begins to travel outward and downward from its point of origin, thereafter encountering and passing through rock units in the subsurface. At each interface between two different rock units, there is the potential for a seismic reflection to occur. The amount of seismic energy that is reflected at an interface is dependent upon the acoustic impedance contrast between the units and the reflection coefficient is one conventional measure of that contrast. The reflection coefficient can be thought of as the ratio of the amplitude of the reflected wave compared with the amplitude of the incident wave. In terms of rock properties: ##EQU1## where, the acoustic impedance of a rock unit is defined to be the mathematical product of the rock density (.rho..sub.1 and .rho..sub.2 being the densities of the upper lower rock units, respectively) multiplied times the velocity in the same rock unit, V.sub.1 and V.sub.2 corresponding to the upper and lower rock unit velocities. (Strictly speaking, this equation is exactly correct only when the wavelet strikes the rock interface at vertical incidence. However, it is generally accepted in the industry that the requirement of verticality is satisfied if the wavelet strikes the interface within about 20.degree. of the vertical.) However, at angles of incidence in excess of about 20.degree., amplitude versus offset effects "AVO", hereinafter) can have a substantial impact on the observed reflectivity, if the reflector is one which might support to an AVO-type reflection.
Reflected energy that is recorded at the surface can be represented conceptually as the convolution of the seismic wavelet with a subsurface reflectivity function: the so-called "convolutional model". In brief, the convolutional model attempts to explain the seismic signal recorded at the surface as the mathematical convolution of the down going source wavelet with a reflectivity function that represents the reflection coefficients at the interfaces between different rock layers in the subsurface. In terms of equations, EQU x(t)=w(t)*e(t)+n(t),
where x(t) is the recorded seismogram, w(t) is the seismic source wavelet, e(t) is the earth's reflectivity function, n(t) is random ambient noise, and "*" represents mathematical convolution. Additionally, the convolutional model requires, in part, that (1) the source wavelet remains invariant as it travels through the subsurface (i.e., that it is stationary and unchanging), and (2) the seismic trace recorded at the surface can be represented as the arithmetic sum of the separate convolutions of the source wavelet with each interface in the subsurface (the principle of "superposition", i.e., that wavelet reflectivity and propagation is a linear system.) Although few truly believe that the convolutional model fully describes the mechanics of wave propagation, the model is sufficiently accurate for most purposes: accurate enough to make the model very useful in practice. The convolutional model is discussed in some detail in Chapter 2.2 of Seismic Data Processing by Ozdogan Yilmaz, Society of Exploration Geophysicists, 1987, the disclosure of which is incorporated herein by reference.
Seismic data that have been properly acquired and processed can provide a wealth of information to the explorationist, one of the individuals within an oil company whose job it is to locate potential drilling sites. For example, a seismic profile gives the explorationist a broad view of the subsurface structure of the rock layers and often reveals important features associated with the entrapment and storage of hydrocarbons such as faults, folds, anticlines, unconformities, and sub-surface salt domes and reefs, among many others. During the computer processing of seismic data, estimates of subsurface velocity are routinely generated and near surface inhomogeneities are detected and displayed. In some cases, seismic data can be used to directly estimate rock porosity, water saturation, and hydrocarbon content. Less obviously, seismic waveform attributes such as phase, peak amplitude, peak-to-trough ratio, and a host of others, can often be empirically correlated with known hydrocarbon occurrences and that correlation a collected over new exploration targets. In brief, seismic data provides some of the best subsurface structural and stratigraphic information that is available, short of drilling a well.
Seismic data are limited, through, in one crucial regard: rock units that are relatively "thin" are often not clearly resolved. In more particular, whereas seismic reflection data can provide a near "geologic cross section" representation of the subsurface when the lithologic layers are relatively "thick", the seismic image that results when the layers are "thin" is much less clear. This phenomenon is known to those skilled in the art as the seismic resolution problem.
Seismic resolution in the present context will be taken to refer to vertical resolution within a single seismic trace, and is loosely defined to be to the minimum separation between two seismic reflectors in the subsurface that can be recognized as separate interfaces--rather than as a single composite reflection--on the seismic record. By way of explanation, a subsurface unit is ideally recognized on a seismic section as a combination of two things: a distinct reflection originating at the top of the unit and a second distinct reflection, possibly of opposite polarity, originating from its base. In the ideal case, both the top and the bottom of the unit appear on the recorded seismogram as distinct and isolated reflectors that can be individually "time picked" (i.e., marked and identified) on the seismic section, the seismic data within the interval between the two time picks containing information about the intervening rock unit. On the other hand, where the seismic unit is not sufficiently thick, the returning reflections from the top and the bottom of the unit overlap, thereby causing interference between the two reflection events and blurring the image of the subsurface. This blurred image of the subsurface is one example of the phenomena known to those skilled in the art as the "thin bed" problem.
FIG. 1 illustrates in a very general way how the thin bed problem arises under the axioms of the convolutional model. Consider first the "thick" bed reflection depicted in FIG. 1A. On the left side of this figure is represented a source wavelet which has been generated on the surface of the earth. The source wavelet travels downward unchanged through the earth along path P1 until it encounters the rock unit interface labeled "A." (Note that the wave paths in this figure are actually vertical, but have been illustrated as angled for purposes of clarity. This is in keeping with the general practice in the industry.) In FIG. 1A, when the down-going seismic waveform encounters Interface "A" a portion of its energy is reflected back toward the surface along path P2 nd is recorded on the surface as the reflected event R1. Note that wavelet R1 is reversed in polarity as compared with the source wavelet, thereby indicating a negative reflection coefficient at the "A" interface. This polarity reversal is offered by way of example only and those skilled in the art recognize that reflection coefficients of either polarity are possible.
The remainder of the down-going energy (after the partial reflection at the interface "A") continues through the thick bed until it strikes Interface "B" at the base of the thick lithographic unit. Upon reaching the "B" interface, part of the wavelet energy continues deeper into the earth along path P5, while the remainder of its energy is reflected back to the surface along path P4 where it is recorded as reflection R2. Note that the reflection from interface "B" occurs at a later point in time than the reflection from interface "A". The exact time separation between the two events depends on the thickness and velocity of the intervening layer between the two interfaces, with thick layers and/or slow velocities creating a greater time separation between the top and base reflections. The temporal thickness of this layer is the time that is required for the seismic waveform to traverse it.
On the surface, the composite thick bed reflection--the expression actually recorded--is the arithmetic sum (superposition) of the two returning reflections, taking into account the time separation between the events. Because the two returning wavelets do not overlap in time, the recorded seismic record clearly displays both events as indicative of two discrete horizons. (Note that the time separation between the two reflected events has not been accurately portrayed in this figure. Those skilled in the art know that the time separation should actually be twice the temporal thickness of the layer.)
Turning now to FIG. 1B, wherein a thin bed reflection is illustrated, once again a source wavelet is generated on the surface of the earth which then travels along path P6 until it encounters the rock unit interface labeled "C." (As before, the wave paths in the figure are actually vertical.) As is illustrated in FIG. 1B, when the down-going seismic waveform encounters Interface "C" a portion of its energy is reflected back toward the surface along path P7, where it is recorded as reflection R3. The remainder of the down-going energy continues through the thin bed until it strikes Interface "D". Upon reaching the "D" interface, part of the wavelet energy continues deeper into the earth along path P10, while the remainder of its energy is reflected back to the surface along path P9 where it is recorded as reflection R4.
Note once again, that the reflection from interface "D" occurs at a later time than the reflection from interface "C", however the temporal separation between the two reflections in the case of a thin bed is less because the distance the waveform must travel before being reflected from interface "D" is less. In fact, the time separation between the two reflections is so small that the returning (upward going) wavelets overlap. Since the composite thin bed reflection is once again the arithmetic sum of the two returning reflections, the actual recorded signal is an event that is not clearly representative of either the reflection from the top or the base of the unit and its interpretation is correspondingly difficult. This indefinite composite reflected event exemplifies the typical thin bed problem.
Needless to say, the thickness of a subsurface exploration target is of considerable economic importance to the oil company explorationist because, other things being equal, the thicker the lithographic unit the greater the 1i it might potentially contain. Given the importance of accurately determining layer thickness, it should come as no surprise that a variety of approaches have been employed in an effort to ameliorate the thin bed problem.
A first technique that is almost universally applied is shortening the length of the seismic wavelet, longer wavelets generally offering worse resolution than shorter ones. During the data processing phase the recorded seismic waveform may often be shortened dramatically by the application of well known signal processing techniques. By way of example, it is well known to those skilled in the art that conventional predictive deconvolution can be used to whiten the spectrum of the wavelet, thereby decreasing its effective length. Similarly, general wavelet processing techniques, including source signature deconvolution and any number of other approaches, might alternatively be applied to attempt to reach a similar end result: a more compact waveform. Although any of these processes could result in dramatic changes to the character of the seismic section and may shorten the length of the wavelet significantly, it is often the case that further steps must be taken.
Even the best signal processing efforts ultimately only postpone the inevitable: no matter how compact the wavelet is, there will be rock layers of economic interest that are too thin for that wavelet to resolve properly. Thus, other broad approaches have been utilized that are aimed more toward analysis of the character of the composite reflection. These approaches are based generally on the observation that, even when there is only a single composite reflection and the thickness of the layer cannot be directly observed, there is still information to be found within the recorded seismic data that may indirectly be used to estimate the actual thickness of the lithographic unit.
By way of example, FIG. 4A illustrates a familiar "pinch out" seismic model, wherein the stratigraphic unit of interest (here with its thickness measured in travel time rather than in length) gradually decreases in thickness until it disappears (i.e., "pinches out") at the left most end of the display. FIG. 4B is a collection of mathematically generated synthetic seismograms computed from this model that illustrate the noise free convolution of a seismic wavelet with the interfaces that bound this layer. Notice that at the right most edge of FIG. 4B, the composite signal recorded on the first trace shows that the reflector is clearly delimited by a negative reflection at the top of the unit and a positive reflection at its base. Moving now to the left within FIG. 4B, the individual reflections at the top and base begin to merge into a single composite reflection and eventually disappear as the thickness of the interval goes to zero. Note, however, that the composite reflection still continues to change in character even after the event has degenerated into a single reflection. Thus, even though there is little direct visual evidence that the reflection arises from two interfaces, the changes the reflections exhibit as the thickness decreases suggest that there is information contained in these reflection that might, in the proper circumstances, provide some information related to the thickness of the thin bed.
The pioneering work of Widess in 1973 (Widess, How thin is a thin bed?, Geophysics, Vol. 38, p. 1176-1180) has given birth to one popular approach to thin bed analysis wherein calibration curves are developed that rely on the peak-to-trough amplitude of the composite reflected thin bed event, together with the peak-to-trough time separation, to provide an estimate of the approximate thickness of the "thin" layer. (See also, Neidell and Poggiagliomi, Stratigraphic Modeling and Interpretation--Geophysical principles and Techniques, in SEISMIC STRATIGRAPHY APPLICATIONS TO HYDROCARBON EXPLORATION, A.A.P.G. Memoir 26, 1977). A necessary step in the calibration process is to establish a "tuning" amplitude for the thin bed event in question, said tuning amplitude occurring at the layer thickness at which maximum constructive interference occurs between the reflections from the top and base of the unit. In theory at least, the tuning thickness depends only on the dominant wavelength of the wavelet, .lambda., and is equal to .lambda./2 where the reflection coefficients on the top and bottom of the unit are the same sign, and .lambda./4 where the reflection coefficients are opposite in sign.
Because of the flexibility of calibration-type approaches, they have been used with some success in rather diverse exploration settings. However, these amplitude and time based calibration methods are heavily dependent on careful seismic processing to establish the correct wavelet phase and to control the relative trace-to-trace seismic trace amplitudes. Those skilled in the seismic processing arts know, however, how difficult it can be to produce a seismic section that truly maintains relative amplitudes throughout. Further, the calibration based method disclosed above is not well suited for examining thin bed responses over a large 3-D survey: the method works best when it can be applied to an isolated reflector on a single seismic line. It is a difficult enough task to develop the calibration curve for a single line: it is much more difficult to find a calibration curve that will work reliably throughout an entire 3-D grid of seismic data. Finally, thin bed reflection effects are occasionally found in conjunction with AVO effects, an occurrence that can greatly complicate thin bed interpretation.
Heretofore, as is well known in the seismic processing and seismic interpretation arts, there has been a need for a method extracting useful thin bed information from conventionally acquired seismic data which does suffer from the above problems. Further, the method should also preferably provide seismic attributes for subsequent seismic stratigraphic and structural analysis. Finally, the method should provide some means of analyzing AVO effects, and especially those effects that occur in conjunction with thin bed reflections. Accordingly, it should now be recognized, as was recognized by the present inventors, that there exists, and has existed for some time, a very real need for a method of seismic data processing that would address and solve the above-described problems.
Before proceeding to a description of the present invention, however, it should be noted and remembered that the description of the invention which follows, together with the accompanying drawings, should not be construed as limiting the invention to the examples (or preferred embodiments) shown and described. This is so because those skilled in the art to which the invention pertains will be able to devise other forms of this invention within the ambit of the appended claims. Finally, although the invention disclosed herein may be illustrated by reference to various aspects of the convolutional model, the methods taught below do not rely on any particular model of the recorded seismic trace and work equally well under broad deviations from the standard convolutional model.